function [x, y, a, blink, speed, vel_x, vel_y, shrink]=EL1k_blinkremover(t, x, y, a, thresh, plot_it)
%[x1 y1 a1 blink speed vel_x vel_y shrink]=EL1k_blinkremover(t, x, y, a, [thresh], [plot_it])
% data must correct to degrees per pixel first.
% inputs:
% t = time in ms
% x = horizontal position in degrees of visual angle
% y = vertical position in degrees of visual angle
% a = area : size of pupil
% thresh = threshold for cutting out 'blink' data and replacing it with interpolated lines
% plot_it = do you want to see my superfluous ploting that takes eons?
%
% ~bcoe:  coe@queensu.ca

%
if nargin==0
    help(mfilename);
    return
end
if nargin<4
    warning ([mfilename ' requires 4 input vectors (t, x, y, a) of equal length']) %#ok<WNTAG>
    help(mfilename);
    return
end
if nargin==4
    thresh=2.5;
    plot_it=0;
end
if nargin==5
    plot_it=0;
end


vel_x=single(zeros(size(x)));
vel_y=single(zeros(size(x)));
speed=single(zeros(size(x)));
shrink=single(zeros(size(x)));
%blink=false(size(x));
if sum(isnan(x))==length(x) || sum(isnan(y))==length(y)
    return
    %             vel_x()=0;
    % 			vel_y()=0;
    % 			speed(starloc(tr):end_loc(tr))=0;
    % 			shrink(starloc(tr):end_loc(tr))=0;
    % 			x(starloc(tr):end_loc(tr))=x;
    % 			y(starloc(tr):end_loc(tr))=y;
    % 			a(starloc(tr):end_loc(tr))=a;
    % 			blink(starloc(tr):end_loc(tr))=1;
    % 			return
end
blink=a<(mean(mean(a(a>0)))/2); % less than half of it's non-zero mean
if isnan(sum(x))
    tdiff=[(t(2)-t(1));(t(3:end)-t(1:end-2))/2;(t(end)-t(end-1))];
    adiff=[(a(2)-a(1));(a(3:end)-a(1:end-2))/2;(a(end)-a(end-1))];
    tshrink=filtfilt([1 1 1]/3,1,double(abs(adiff./tdiff))); % filtfilt now requires doubles
    tshrinkstd=std(tshrink(tshrink<10))*thresh*2;% temp var just to shore up the ends
    firstgood=find(tshrink<tshrinkstd  &~blink & ~isnan(x) ,1,'first');
    lastgood=find(tshrink<tshrinkstd  &~blink & ~isnan(x) ,1,'last');
    if firstgood>1
        x(001:firstgood)=x(firstgood);
        y(001:firstgood)=y(firstgood);
        a(001:firstgood)=0;
    end
    if lastgood<length(x)
        x(lastgood:end)=x(lastgood);
        y(lastgood:end)=y(lastgood);
        a(lastgood:end)=0;
    end
    %a(1)=a(find(isnan(a)==0,1,'first'));
    %a(end)=a(find(isnan(a)==0,1,'last'));
end %isNaN
blink=a<(mean(mean(a(a>0)))/2); % less than half it's non-zero mean
[speed, vel_x, vel_y, shrink]=xyt2speed(x,y,t,a);
gooddata=~isnan(speed);
speed(gooddata)=single(filtfilt([1 1 1]/3,1,double(speed(gooddata))));% filtfilt now requires doubles
gooddata=~isnan(shrink);
shrink(gooddata)=single(filtfilt([1 1 1 1]/4,1,double(shrink(gooddata))));% filtfilt now requires doubles

if plot_it
    ttt=t-t(1); %just for plotting
    figure(1);clf; hold on;
    ylim([-40 20])
    xlim([0 ttt(end)])
    plot(ttt,abs(tshrink)-30,'color',[1 .5 0],'linewidth',1.5);
        plot(ttt,a/50-30,'k','linewidth',2)
    plot(ttt,speed/50-30,'g','linewidth',2);
plot(ttt,x,'b','linewidth',2)
    plot(ttt,y,'r','linewidth',2)
    %hline(100/50-30,'k'); % 50 dps
    %hline(100/50-30,'--g'); % 50 dps
    hline((mean(mean(a(a>0)))/2)/50-30,'--k'); % dynamic thres
end

s_blink01=filtfilt([1 1 1 1 1 1 1]/7,1,double(blink));% filtfilt now requires doubles
s_blink02=filtfilt([1 1 1 1 1 1 1]/7,1,round(s_blink01+.25));% filtfilt now requires doubles
blink=logical(round(s_blink02+.25));
blink(end)=0;
%loc=find(diff([0; blink]))-1;
loc=(diff([0; blink]));
loce=[find(loc==1) find(loc==-1)-1];

spinner={'\','|','/','-'};
fprintf('running blink remover:-\n');
if ~isempty(loc) %actually blinked
    shrinkstd=std(shrink(shrink<10))*thresh;
    speedstd=nanstd(speed(speed<200 & speed>0))*thresh;
    if plot_it
        plot(ttt(blink),a(blink)/50-30,'xk');
        vfill(ttt(loce(:,1)),ttt(loce(:,2)),[],.2 )
        %vfill(ttt(loc(1:2:end)),ttt(loc(2:2:end)),[],.2 )
        vline(ttt(loce))
        hline(shrinkstd-30,'c');
        hline(shrinkstd-30,[1 .5 0]);
        hline(speedstd/50-30,'-.g');
    end
    
    sloc=[1; loce(:,1);length(blink)];
    eloc=[1; loce(:,2);length(blink)];
    for i = 2:length(sloc)-1
        fprintf('\b\b%s\n',spinner{mod(i,4)+1});
        %forward in time
        eee=find(speed(eloc(i):sloc(i+1))<speedstd &...
            abs(shrink(eloc(i):sloc(i+1)))<shrinkstd ...
            , 1)+eloc(i)-1 ;
        if isempty(eee)% try again with higher threshold
            eee=find(speed(eloc(i):sloc(i+1))<speedstd*1 &... % ok *1 means do nothing
                abs(shrink(eloc(i):sloc(i+1)))<shrinkstd*1 ...
                , 1)+eloc(i)-1 ;
            if isempty(eee)% % try with next saccade
                eee=sloc(i+1);
                sloc(i+1)=sloc(i);
            end
        end
        blink(eloc(i):eee)=1;
        if plot_it && eee
            xlim([ttt(eloc(i))-50 ttt(eloc(i))+150 ])
            vfill(ttt(eloc(i)), ttt(eee),[],.5)
            %drawnow
        end
        eloc(i)=eee;% redefine end
        
        % backward in time
        sss=find(speed(eloc(i-1):sloc(i))<speedstd & ...
            abs(shrink(eloc(i-1):sloc(i)))<shrinkstd ...
            ,1,'last')+eloc(i-1)-1;
        if isempty(sss)% try again with higher threshold
            sss=find(speed(eloc(i-1):sloc(i))<speedstd*1 & ...% ok *1 means do nothing
                abs(shrink(eloc(i-1):sloc(i)))<shrinkstd*1 ...
                ,1,'last')+eloc(i-1)-1;
            if isempty(sss)% try with next saccade
                sss=sloc(i-1);
                %sloc(i+1)=sloc(i-1);
            end
        end
        if sss-eloc(i-1)<5 && sss~=sloc(i-1); % if gap of real data too small
            drawnow
            sss=sloc(i-1);
        end
        blink(sss:sloc(i))=1;
        if plot_it
            xlim([ttt(sloc(i))-150 ttt(sloc(i))+50 ])
            vfill(ttt(sss), ttt(sloc(i)),[],.5)
            %drawnow
        end
        sloc(i)=sss;% redefine start

        % interpolate
        %try
        if length(sss:eee)>500
            drawnow
        end
        x(sss:eee)=linspace(x(sss),x(eee),eee-sss+1);
        y(sss:eee)=linspace(y(sss),y(eee),eee-sss+1);
        %a(sss:eee)=linspace(a(sss),a(eee),eee-sss+1);
        % recalculate
        [speed, vel_x, vel_y, shrink]=xyt2speed(x,y,t,a);
        %speed=single(filtfilt([1 1 1]/3,1,double(speed)));% filtfilt now requires doubles
        %shrink=single(filtfilt([1 1 1 1]/4,1,double(shrink)));% filtfilt now requires doubles
        gooddata=~isnan(speed);
        speed(gooddata)=single(filtfilt([1 1 1]/3,1,double(speed(gooddata))));% filtfilt now requires doubles
        gooddata=~isnan(shrink);
        shrink(gooddata)=single(filtfilt([1 1 1 1]/4,1,double(shrink(gooddata))));% filtfilt now requires doubles
        
        if plot_it
            % replot
            xlim([ttt(sss)-150 ttt(eee)+150 ] )
            plot(ttt(sss:eee), x(sss:eee),'ob')
            plot(ttt(sss:eee), y(sss:eee),'or')
            plot(ttt(sss:eee), a(sss:eee)/50-30,'ok')
            plot(ttt(sss:eee), speed(sss:eee)/50-30,'og')
            plot(ttt(sss:eee), abs(shrink(sss:eee))-30,'o','color',[1 .5 0])
            drawnow
            pause (.1)
        end
        %catch %#ok<CTCH>
        %	drawnow;
        %end
    end % sloc
end % actually blinked

end % function



